PsL Monthly 1993 December/PSL Monthly Shareware CD-ROM (December 1993).iso/prgmming/dos/c/snippet.exe/RAYS.C

C/C++ Source or Header | 1991-09-23 | 5.0 KB | 180 lines

/*
**  public domain ray tracing functions by Daniel Lyke
*/

#include <math.h>
#include "plot.h"


struct RAY
{
      double dx, dy, dz; /* Direction vector */
      double ox, oy, oz; /* Origin */
};

struct PLANE
{
      double nx, ny, nz; /* Vector normal (perpendicular) to plane */
      double px, py, pz; /* Point on plane */
};

struct SPHERE
{
      double cx, cy, cz; /* Center of sphere */
      double r2; /* Radius squared */
};

struct VECTOR
{
      double dx, dy, dz; /* Three dimensional vector */
};


/*
** Return the closest point of intersection of a ray with a sphere.
** Rays starting within the sphere are returned as not intersecting.
*/

double sphere_intersect(struct RAY ray, struct SPHERE sphere)
{
      double a, b, c, t1, t2, t3, close, farther;

      a = ray.dx * ray.dx + ray.dy * ray.dy + ray.dz * ray.dz;

      close = farther = -1.0;

      if(a)
      {
            b = 2.0 * ((ray.ox - sphere.cx) * ray.dx
                  + (ray.oy - sphere.cy) * ray.dy
                  + (ray.oz - sphere.cz) * ray.dz);
            c = (ray.ox - sphere.cx) * (ray.ox - sphere.cx)
                  + (ray.oy - sphere.cy) * (ray.oy - sphere.cy)
                  + (ray.oz - sphere.cz) * (ray.oz - sphere.cz) - sphere.r2;

            t1 = b * b - 4.0 * a * c;
            if(t1 > 0)
            {
                  t2 = sqrt(t1);
                  t3 = 2.0 * a;
                  close = -(b + t2) / t3;
                  farther = -(b - t2) / t3;
            }
      }
      return (double)((close < farther) ? close : farther);
}


/*
** Return time to intersection with a plane
*/

double plane_intersect(struct RAY ray, struct PLANE plane)
{
      double p1, p2, p3;

      p1 = plane.px * plane.ny + plane.py * plane.ny + plane.pz * plane.nz;
      p2 = ray.ox * plane.nx + ray.oy * plane.ny + ray.oz * plane.nz;
      p3 = ray.dx * plane.nx + ray.dy * plane.ny + ray.dz * plane.nz;

      return (double)((p1-p2)/p3);
}


/*
** Given a ray normal to the surface of reflection (n), an incident ray (i),
** return a reflected ray (r)
*/

void reflect(struct VECTOR *n, struct VECTOR *i, struct VECTOR *r)
{
      double ndotn, idotn;

      ndotn = (n->dx * n->dx + n->dy * n->dy + n->dz * n->dz);
      idotn = (n->dx * i->dx + n->dy * i->dy + n->dz * i->dz);

      r->dx = i->dx - (2.0 * (idotn) / ndotn) * n->dx;
      r->dy = i->dy - (2.0 * (idotn) / ndotn) * n->dy;
      r->dz = i->dz - (2.0 * (idotn) / ndotn) * n->dz;
}

int trace(void)
{
      int x,y,c;
      static struct PLANE plane = { 0.0, 1.0, 0.001, -6.0, 0.0, 0.0};
      static struct SPHERE sphere = { 0.0, 2.0, 8.0, 49.0 };
      struct RAY ray;

      struct VECTOR v1, v2, v3;

      double t1, t2, time;

      c=0;

      for(x=0;x<320 && !kbhit(); x+=2)
      {
            for(y = 0; y < 200; y+=2)
            {
                  /* Ray needs to be reset because of
                     alterations in reflections       */

                  ray.ox = 0.0;
                  ray.oy = 0.0;
                  ray.oz = 0.0;

                  ray.dz = 1.0;
                  ray.dy = -((double)y - 100.0) / 75.0;
                  ray.dx = ((double)x - 160.0) / 80.0;

                  t1 = sphere_intersect(ray,sphere);
                  t2 = plane_intersect(ray,plane);

                  if(t1 > 0.0 && (t2 < 0.0 || t2 > t1)) /* Circle in fore? */
                  {
                        v1.dx = ray.dx; v1.dy = ray.dy; v1.dz = ray.dz;

                        v2.dx = ((ray.dx * t1 + ray.ox) - sphere.cx);
                        v2.dy = ((ray.dy * t1 + ray.oy) - sphere.cy);
                        v2.dz = ((ray.dz * t1 + ray.oz) - sphere.cz);

                        reflect(&v2,&v1, &v3);
                        ray.ox += ray.dx * t1;
                        ray.oy += ray.dy * t1;
                        ray.oz += ray.dz * t1;
                        ray.dx = v3.dx;
                        ray.dy = v3.dy;
                        ray.dz = v3.dz;
                        t2 = plane_intersect(ray,plane);
                        if(t2 > 0.0)
                        {
                              int cr;

                              cr = (abs((int)(t2 * ray.dz + ray.oz)) % 2) +
                                    2 * (abs((int)(t2 * ray.dx + ray.ox))
                                    % 2);
                              plot(x,y,cr);
                        }
                        else  plot(x,y,1);
                  }
                  else if(t2 > 0.0)
                  {
                        int cr;

                        cr = (abs((int)(t2 * ray.dz + ray.oz)) % 2) + 2 *
                              (abs((int)(t2 * ray.dx + ray.ox)) % 2);
                        plot(x,y,cr);
                  }
            }
      }
}

void main(void)
{
      set_mode(4); /* Built for CGA, 4 color mode through BIOS */

      trace();

      getch();
      set_mode(3);
}